Physics, Stability and Dynamics of Supply Networks 
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We show how to treat supply networks as physical transport problems governed by balance equa- 
tions and equations for the adaptation of production speeds. Although the non-linear behaviour is 
different, the linearized set of coupled differential equations is formally related to those of mechan- 
ical or electrical oscillator networks. Supply networks possess interesting new features due to their 
complex topology and directed links. We derive analytical conditions for absolute and convective 
instabilities. The empirically observed "bull-whip effect" in supply chains is explained as a form 
of convective instability based on resonance effects. Moreover, it is generalized to arbitrary supply 
networks. Their related eigenvalues are usually complex, depending on the network structure (even 
without loops). Therefore, their generic behavior is characterized by oscillations. We also show that 
regular distribution networks possess two negative eigenvalues only, but perturbations generate a 
spectrum of complex eigenvalues. 
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Econophysics has stimulated a lot of interesting re- 
search on problems in finance and economic systems Q, 
applying methods from statistical physics and the the- 
ory of complex systems. Recently, the dynamics of sup- 
ply networks 0-LjHil]- he. the flow of materials through 
networks, has been identified as an interesting physical 
transport problem [J- |J- LJ. which is also reflected by 
the notion of "factory physics" [|| . Potential applications 
reach from production networks to business cycles, from 
metabolic networks to food webs, up to logistic prob- 
lems in disaster management. Empirical and theoretical 
studies have shown that the flow of goods between differ- 
ent producers or suppliers can be described analogously 
to driven many-particle flows between sources and sinks 
(depots), where the particles represent materials, goods, 
or other resources. In contrast to the stationary behav- 
ior assumed by most production engineers, this flow may 
display a complex dynamics in time, including oscillatory 
patterns and chaos pi |8j • A particular focus has been on 
the empirically observed and well-known "bull-whip ef- 
fect" [il |H IS lU ; which describes the amplification 
of the oscillation amplitudes of delivery rates in supply 
chains compared to the variations in the consumption 
rate of goods. 

A promising approach to the non-linear interactions 
and dynamics of supply chains is based on "fluid-dynamic 
models" OaM, which are related to macroscopic traffic 
models jj, 0, |5j In contrast to classical approaches like 
queueing theory and event-driven simulations, they are 
better suited for an on-line control under dynamically 
changing conditions. These fluid-dynamic models have 



recently been generalized to cope with discrete spaces 
(production steps) 0, 0|. However, the impact of the 
topology of supply networks, connecting the subject to 
the statistical physics of networks 0] nas n °t Y e t been 
thoroughly investigated. Its relevance for the stability 
and dynamics of supply networks will, therefore, be the 
main focus of this paper. Here, we will specifically con- 
centrate on analytical results for the dynamic behavior 
in the linear regime around the usually assumed station- 
ary state, while non-linear effects are investigated else- 
where The corresponding equations can be mapped 
onto the ones of particular mechanical or electrical net- 
works, but supply networks are generally more complex 
and their non-linear dynamics is different. In particular, 
supply networks are directed and possess other charac- 
teristic topologies 0]. Hence, our analytical investiga- 
tion yields interesting new results. Apart from a gener- 
alization of the bull-whip effect from linear (sequential) 
supply chains to networks, they allow us to explain vari- 
ous surprising features of numerical simulations of supply 
networks: (i) The bull-whip effect can occur even though 
the eigenvalues of linear supply chains are negative, (ii) 
Supply networks tend to show oscillations, even without 
loops in the material flows, (iii) Regular distribution net- 
works are characterized by two negative eigenvalues, but 
random perturbations in the network structure cause a 
spectrum of oscillating modes. 

Our simplified model of supply networks consists of 
u suppliers i delivering products to other suppliers j or 
consumers. The suppliers j (e.g. production units of a 
plant or company) will be assumed to deliver products, 
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resources or materials k at a certain rate d k jXj(t). Their 
demand of goods of kind i from other suppliers j is given 
by CijXj(t), where < cy < 1. The proportionality 
CijXj(t) to the production rate Xj(t) of products j re- 
flects that a quantity Cy of product i is needed to produce 
one unit of product j. The coefficients Cij define an in- 
put matrix C = (cij) and dij an output matrix D. The 
interactions among the different suppliers i and j are rep- 
resented by the supply matrix S = D — C. In addition, 
we need to consider the inventories, stock levels, or avail- 
ability Ni(t) of products, materials, or resources i. 

Our model for the dynamics of supply networks con- 
sists of two sets of equations, one for the change of in- 
ventories N (t) with time t and another one for the adap- 
tation of the delivery or production rates Xj(t) of the 
suppliers or production units j. The inventories Ni(t) of 
goods of kind i are described by material balance equa- 
tions for the flows of goods, which could be viewed as a 
discontinuous version of the continuity equation: 



inflow 



re — entrant 
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dt 



J2 d lk x k (t) - J2 cy-Xj-C*) + Y i(t) 

k=i 



(1) 
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with < Cij^dij < 1 (1 < i,j < u) and the "normal- 
ization conditions" dio = 1 — J2t=i — 0, (H,u+i = 
1 — Y]j—i Cij > 0. Frequently, one assumes dij — Sij with 
Sij = 1 for i = j and 5ij = otherwise [||. This corre- 
sponds to situations in which production units or suppli- 
ers are defined such that they deliver one kind of good 
only. The quantity 

Yi{t)= Ci, u+1 X +1 (t) - (ImXoW (2) 

consumption and losses inflow of resources 

comprises the consumption rate of goods i, losses, and 
waste (the "export" of material) , minus the inflows into 
the considered system (the "imports"). 

Changes in the demand Yi(t) sooner or later require 
an adaptation of the production rates Xj (t) . This adap- 
tation takes time (to determine the change in demand, 
to finish the planning process and administrative steps, 
and to adapt orders and production capacities). We will 
reflect this by an adaptation time T and assume that 
the change in the production rate is proportional to 
the deviation of the actual production rate Xj from 
the desired one Wj({N}, {dNi/dt}), which depends on 
the inventories N and their changes dNi/dt in time: 
dXj/dt = [Wj({Ni}, {dNi/dt}) - Xj{t)]/T. Although a 
more general treatment is possible, we will first focus on 
cases where the delivery rate Xj of supplier j is only 
adapted to the inventory N of the dominant product 
i given by dij = max^ d^j ■ Moreover, we will enumer- 
ate suppliers j according to their dominating products i, 



which implies j = i. The resulting adaptation equation 

^ = i [Wi (N, dNi/dt) - Xi(t)} (3) 

still implies an indirect dependence on the production 
and delivery rates Xj of other suppliers j ^ i via the 
dependence on dN/dt, see Eq. Q). 

For the desired production rate Wi(N, dN /dt),we will 
assume that (i) it is non-negative and decreasing with in- 
creasing inventories, (ii) it is proportional to the steady- 
state inventory Ni, (iii) it responds to the relative devi- 
ation Ni/N of the inventory Ni from the steady-state 
inventory N and the relative change (dN / dt) / N(t) in 
time. This suggests the scaling relation 



Wi N 



dN 
dt 



— {Ni 1 dN, 

vNiP =, 

\N N dt 



( 4 ) 



with a scaling function P which, for simplicity, is assumed 
to be identical for all sectors. The stationary value N of 
the inventory Ni can be easily determined by solving the 
linear system of equations 

u u 
Yi = ^(dy - Cij )Xj =VP(1, 0) £)(dy - Cij )Nj , (5) 



where Yi denotes the time-average of the demand Yi(t), 
while Xi and N t denote the stationary values correspond- 
ing to the case Yi — Yi. By appropriate choice of v, we 
can set P(1,0) = 1. 

Since we want to know whether the stationary state 
of the supply network is stable, as is often implicitly as- 
sumed in queueing theory, we will focus on what hap- 
pens in case of small deviations Xj(r) = [Xi(rT) — X{\T, 
ni(r) = Ni(rT)-Ni, and w (t) = [^(tT) -Y t ]T from it. 
Here, we have introduced a scaling of the model variables 
to dimensionless units. For example, r = t/T is a dimcn- 
sionless time, i.e. in the following we will measure the 
time in units of T. Close to equilibrium, we can linearize 
the model equations. In matrix notation they read 



d_ 

(It 



u(t) 
x{t) 



= M 



with M = 



n(r) 
x(t) 

, 
AE , 



-v(t) 
By(r) 

S 

E-BS 



(6) 
(7) 



Here, A = -vTdP(l, 0)/dzi and BT 
—i/TdP(l,0)/dz2 denote the negative values of the 
partial derivatives of TWi(Ni, dNi/dt) with respect 
to the first variable Z\ — N/N and the second one 
z 2 = (dN/dt) /Ni = (dNi/dr)/ (NiT) in the equilibrium 
point (N, dN/dt) = (N,0). The negative signs reflect 
that A and B are typically positive, as the production 
rate should decrease when the inventory increases. 

The above linear system of coupled differential equa- 
tions is a quite general approach to the study of supply 
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networks and material flows close to the stationary state. 
It describes the response of delivery rates to a variation 
y(t) in the consumption rate and can be derived from 
various non-linear supply chain models which differ in 
their degree of detailedness concerning the consideration 
of forecasts 0, price dynamics 0,0]) or lack of materials 
Therefore, it is worth investigating the dependence 
on the dimensionless parameters A and B, and on the 
structure of the supply network characterized by S. The 
system of 2u differential equations can also be written as 
a system of u second-order differential equations 



+ [E + BS]^ + ASx(r)=5(t) 

> dy 



(8) 



where g(t) — Ay{r) + B^. Note that there exists a 
matrix T which allows one to transform the matrix 
E - S via T -1 (E - S)T = J into either a diagonal or 
a Jordan normal form J. Defining p(r) — T _1 x(r) and 
h(r) = T _1 <?(t), we obtain the coupled set of second- 
order differential equations 



d Pi dm 2 

~ j q - + 27i-^ h u>i m 

ar* dr 



dp 



h [Bp l+1 + A^-^j + hi , 

(9) 

where 7i = [1 + B(l - J K )]/2, w< = [A(l ~ J,*)] 1/2 , and 
h = Ji.i+i- This can be interpreted as a set of equa- 
tions for linearly coupled damped oscillators with damp- 
ing constants 74, eigenfrequencies u>i, and external forcing 
hi(r). The other forcing terms on the right-hand side are 
due to interactions of suppliers. They appear only if J is 
not of diagonal, but of Jordan normal form. Because of 
bi = Ji.i+i, Eqs. © can always be analytically solved in 
a recursive way, starting with the highest index i = u. 
Note that, in the case D = E (i.e. dij — Sij), J a are 
the eigenvalues of the input matrix C and < \Ju\ < 1. 

Equation JjJJl has a special periodic solution of the form 
= M o e i(at-ft) j ^(f) = h o e i a t^ where j = de _ 

notes the imaginary unit. Inserting this into JJjJ and di- 
viding by e lat immediately gives 

(-a 2 +2ia-/ l +uj l 2 )fi a l e^ = b i (A+iaB)p° i+1 e- i0i + l +h° i . 

(10) 

With e^ = cos(0) ± isin(0) this implies 



y/Re 2 + Im 2 &' 



x /(wi 2 -a 2 ) 2 + (2a 7i ) 2 e 1 ^ ' 



(11) 



where Re = bip° +1 [Acos((3i + i) + aB sin(/?j+i)] + h® and 
Im = bii4 +1 [aB cos(/3 i+ i) - Asm(/3 i+1 )]. 



l lA 2 + (aB) 2 }M +1 ) 2 + h°H l + (h°y 

Mi ~ V - « 2 ) 2 + (2«7,) 2 1 ' 

with Hi = 2bi/i° +1 [Acos(P l+ i) + aB sm((3 i+ i)]. Finally, 
we have (3i — ifi — pi with tan^ = 2aji/(uji 2 — a 2 ) and 

birf +1 [aB cos(/? l+ i) - Asin(A+i)] 



tan pi 



bipP i+1 [Acos(/3 i+1 ) + aBsm(f} l+1 )} + h° 



(13) 



For /i° = 0, we obtain tan(< / 9 i — /?,) = tan(5 — Pi+i) with 
tan S = aB/A, i.e. the phase shift between i and i + 1 is 
just (3i - = ifi - 8. 

According to Eq. , the dynamics of our supply net- 
work model can surprisingly be reduced to the dynamics 
of a linear (sequential) supply chain, but with the new 
entities i having the meaning of "quasi-suppliers" (anal- 
ogously to "quasi-species" |lj3, |l4j) defined by the lin- 
ear combination p(r) = T x(t). This transformation 
makes it possible to define the bull-whip effect for arbi- 
trary supply networks: It occurs if the amplitude p® is 
greater than the oscillation amplitude p® +1 of the next 
downstream supplier i+1 and greater than the amplitude 
h® of the external forcing, i.e. if p®/ max(p® +1 , h®) > 1. 
This resonance effect corresponds to the case of convec- 
tive instability. According to formula l|12|) , the bull- whip 
effect is particularly likely to appear, if the oscillation 
frequency a of the consumption rate is close to one of 
the resonance frequencies 10 i. 




FIG. 1: (a) Empirical example for a temporal amplification 
in the oscillation amplitude of the weekly order flow of a ma- 
jor company around its average (7 = 0.05, to = 2%/T with 
T = 2 weeks) . (b) Example of a supply network (regular dis- 
tribution network) with two degenerated real eigenvalues (see 
squares in the right subfigure). (c) The eigenvalues A,:,± of 
the randomly perturbed supply network (crosses) are mostly 
complex and located within Gersgorin's discs (large circles). 

Depending on the respective network structure, it can 
also happen that the oscillation amplitude of p{t) is am- 
plified in the course of time (see Fig. la). This case of 
absolute instability can occur if at least one of the eigen- 
values Xi.± of the homogeneous equation © resulting 
for hi — bi — has a positive real part. The (up to) 2u 
eigenvalues 



A 



z,± 



-7i ± vV^i 



(14) 



depend on the (quasi-) supplier i and determine the tem- 
poral evolution of the amplitude of deviations from the 
stationary solution. One can distinguish various interest- 
ing cases: (i) If the supply matrix S is symmetric, as 
for most mechanical or electrical oscillator networks, all 
eigenvalues Ju are real. Consequently, if uji < 7; (i.e. if 
A is small enough), the eigenvalues Aj,± of M are real 
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and negative, corresponding to an overdamped behavior. 
However, if B is too small, the system behavior may be 
characterized by damped oscillations, (ii) Most natural 
and man-made supply networks have directed links, and 
S is not symmetric. Therefore, some of the eigenvalues 
Ju will normally be complex, and an overdamped behav- 
ior is untypical. The characteristic behavior is rather of 
oscillatory nature (although asymmetry does not always 
imply complex eigenvalues |12j). For small values of B 1 
it can even happen that the real part of an eigenvalue 
\ it ± becomes positive. This implies an amplification of 
the oscillations in time (until the oscillation amplitude is 
limited by nonlinear terms) . Surprisingly, this also applies 
to most upper triangular matrices, i.e. when no loops in 
the material flows exist, (iii) Another relevant case are 
linear (sequential) supply chains and regular supply net- 
works (see, e.g., Fig. lb). These are mostly character- 
ized by degenerated zero eigenvalues Ju = and Jordan 
normal forms J, i.e. non-vanishing upper-diagonal ele- 
ments Ji t i+i- Sequential supply chains and regular dis- 
tribution systems belong to this case, which is charac- 
terized by the two it-fold degenerate eigenvalues X± = 
-(l+B)/2± + B) 2 /A - A independently of the sup- 
pliers i. For small enough values A < (1 + £?) 2 /4, these 
systems show overdamped behavior, otherwise damped 
oscillations. 

Note that already very small perturbations in the net- 
work structure can qualitatively change the dynamics of 
supply networks. In case (iii), for example, they cause 
a cluster of complex eigenvalues around the two neg- 
ative eigenvalues A± (if A is small, see Fig. lc). How 
can we explain this interesting observation? We may use 
Gersgorin's theorem on the location of eigenvalues |l5| . 
Applying it to the perturbed matrix B,, — E — S + rjP, 
for small enough values of r), the corresponding eigenval- 
ues should be located within discs of radius rjj^j \Pi^\ 
around the eigenvalues Ju of the original matrix E — S. 
The parameter rj with < r\ < 1 allows to control the 
size of the perturbation. Moreover, P( r ') = Rl^PR^, 
where R,, is the matrix which transforms B,, to a di- 
agonal matrix T>^\ i.e. R^B^R,, = DW. (This as- 
sumes a perturbed matrix B,, with no degenerate eigen- 
values.) Similar discs as for the eigenvalues of B r) can be 
determined for the associated eigenvalues ± of the per- 
turbed (2u x 2u)-matrix M n belonging to the perturbed 
u x u-matrix B^, see Eq. (J7J) and Fig. lc. 

In summary, this contribution has shown how supply 
systems can be modeled by equations for material flows 
in networks. We could explain the empirically observed 
bull-whip effect as a convective instability phenomenon 
based on a resonance effect and generalize it from lin- 
ear (sequential) supply chains to arbitrary supply net- 
works, as these can be transformed on supply chains for 
"quasi-suppliers" . However, the eigenvalues may become 
complex. The dynamics of supply networks depends very 
sensitively on their topology and two parameters A and 



B (which are related to derivatives of the "control func- 
tion" of the delivery rates). In contrast to most systems of 
mechanical or electrical oscillators, supply networks have 
directed links and are asymmetric. Therefore, most sup- 
ply networks are expected to show an damped or growing 
oscillatory behavior, even without loops. Although neg- 
ative eigenvalues are found for regular distribution net- 
works and supply chains, already small perturbations of 
the network structure will cause a spectrum of complex 
eigenvalues, i.e. a qualitative change in the dynamics. Our 
current research focusses on the application to empirical 
input-output data and business cycles |l2| , on non-linear 
properties of supply networks, and applications to disas- 
ter management with a lack of resources. 
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